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Single-molecule force spectroscopy experiments, as well as a number of other physical systems, are 
governed by thermally activated transitions out of a metastable state under the action of a steadily 
increasing external force. The main observable in such experiments is the distribution of the forces, 
at which the escape events occur. The challenge in interpreting the experimental data is to relate 
them to the microscopic system properties. We work out a maximum likelihood approach and show 
that it is the optimal method to tackle this problem. When fitting actual experimental data it is 
unavoidable to assume some functional form for the force-dependent escape rate. We consider a 
quite general and common such functional form and demonstrate by means of data from a realistic 
computer experiment that the maximum number of fit parameters that can be determined reliably 
is three. They are related to the force- free escape rate and the position and height of the activation 
barrier. Furthermore, the results for the first two of these fit parameters show little dependence on 
the assumption about the manner in which the barrier decreases with the applied force, while the 
last one, the barrier height in the absence of force, depends strongly on this assumption. 

PACS numbers: 82.37.Np, 33.15.Fm, 02.50.-r 



I. INTRODUCTION 

A quite remarkable experimental achievement of the 
last decade is the direct observation of chemical disso- 
ciation at the single-molecule level by applying time- 
dependent external forces on the pico-Newton scale. This 
technique is called dynamic force spectroscopy or single- 
molecule force spectroscopy and reviewed e.g. in [HQ- 
It allows to extract kinetic constants and energy land- 
scape parameters of various interactions like antibody- 
antigen recognition Q or protein-DNA interactions 
Also the dynamics of various other experimental systems 
are governed by thermally activated transitions out of a 
metastable state over a potential barrier, whose height 
decreases in time due to a steadily increasing exter- 
nal force. Examples include the polarization reversal in 
nanomagnets friction at the atomic scale and 
phase jumps in Josephson junctions ^, H^. In all these 
cases, the knowledge of the force-dependent rate out of 
the metastable state can be exploited to characterize the 
system studied. The main theme of the present work is 
how to perform such a characterization in the most op- 
timal way and to point out the limitations even of such 
an optimized procedure under practical conditions. 

A typical experimental setup is schematically sketched 
in Fig. 1: the single chemical bond of interest, e.g. in 
a ligand-receptor complex, is connected via two linker 
molecules with the tip of an AFM (atomic force micro- 
scope) cantilever (or some other micromechanical tool) 
and a piezoelectric element. The latter is employed for 
"pulling down" the attached linker molecule at some 
constant velocity, leading to an elastic reaction force of 
the cantilever, determined from the deflection of a laser 
beam. The main quantity of interest is the magnitude of 
the force at the moment when the bond breaks. 

The theoretical interpretation of the observed rupture 
forces is a non-trivial task for the following reasons. Upon 
repeating the same experiment with the same pulling ve- 




FIG. 1: (Color online) Schematic illustration of a single- 
molecule force spectroscopy experiment. Receptor and ligand 
are connected via suitable linkers to the surface and the AFM 
tip, respectively. The distance of the tip form the surface 
can be controlled with a piezoelectric element (not shown). 
When pulled down at constant speed an (approximately lin- 
early) increasing force acts on the bond which can be mea- 
sured through the deflection of the cantilever. 



locity, the rupture forces are found to be distributed over 
a wide range, contrary to what one would naively ex- 
pect for a purely mechanical breaking of a compound 
object at some fixed, "critical" strain force. A further 
theoretical challenge represents the experimental find- 
ing that for different pulling velocities different such dis- 
tributions are obtained. Hence, neither a single rup- 
ture event nor the average rupture force at any fixed 
pulling velocity can serve as a meaningful characteris- 
tics of a given chemical bond strength. Major steps in 
solving the puzzle are due to Bell jTl| and to Evans 
and Ritchie recognizing that a forced bond rup- 
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FIG. 2: Schematic illustration of the total relevant potential 
energy landscape U{x) — xf of a receptor- ligand-bond as func- 
tion of the reaction coordinate x without and with externally 
applied bias forces /i < /2. For low forces the main effect 
is that the energy barrier Eh{f) is lowered by an amount 
Ai?b — Xbf, where Xh is the distance between potential well 
and barrier at zero force. For larger forces this distance de- 
creases. 

tiire event is a thermally activated decay of a metastable 
state that can be described within the general frame- 
work of Kramers reaction rate theory . Subsequently, 
their basic theoretical approach has been extended and 
refined in several important directions, see e.g . Refs. 
d, 2, 14, 15, 16, 17, 18, 19, 20, 2L 22^ 23i,[2l[lf SIH. 

Following Evans and Ritchie [ij, |T2| , a single-molecule 
dissociation process is viewed as thermally activated es- 
cape of a reaction coordinate x over a potential barrier, 
see Fig. 2. Given the one-dimensional potential land- 
scape along the reaction coordinate, the dissociation rate 
k{f) for an instantaneous force / (projected onto the di- 
rection of the reaction coordinate) can be written accord- 
ing to Kramers reaction-rate theory in the form 

k{f)^iu{f)exp{-~Eb{f)/kBT), (1) 

where the pre-exponential factor has the intuitive 

meaning of an "attempt frequency" and the exponen- 
tially leading Boltzmann-Arrhenius factor contains the 
relevant potential barrier £'(,(/) against escape, Boltz- 
mann's constant ks, and the temperature T. 

In dynamic force spectroscopy experiments, the rate 
at which the force increases is much slower than all other 
relevant molecular relaxation processes. Due to this sep- 
aration of time scales , the reaction kinetics can be 
very accurately approximated by the following first-order 
differential equation 

^ = -k{m)n{t), (2) 

where n{t) denotes the survival probability of the bond 
up to time t and k{J{t)) is the accompanying dissociation 
rate ([1]) at an instantaneous external force f{t). 



Ideally, one would wish to use the experimentally es- 
tablished force-dependent escape rate fc(/) to reconstruct 
the potential landscape, i.e. potential energy vs. the re- 
action coordinate, cf. Fig. [51 However, in view of the 
fact that the escape rate depends on the energy differ- 
ence at two force-dependent extrema, this problem does 
not have a unique solution '28*1 . Therefore, one may start 
with some model energy landscape, and try to deduce its 
global features, such as barrier height in the absence of 
the force and dissociation length (the distance between 
potential well and barrier without bias force, cf. Fig. [2]). 
This introduces some specific functional form for the es- 
cape rate ([T]) involving several parameters, which are 
then determined by fitting the experimental data. In the 
present work, we describe the application of the maxi- 
mum likelihood approach [2^, [s^l as a tool to deduce the 
model parameters. We show that this method is superior 
to any other approach that may be used for this purpose. 

After introducing the method, we discuss its applica- 
tion to determine the parameters of a commonly used 
Ansatz for the rate due to Bell, henceforth called Bell's 
model [ill, This approximation assumes that the 

force is sufficiently small, so that the force-dependent 
barrier height decreases linearly with the force, the pro- 
portionality constant Xh being the dissociation length 
(see Fig.H]), 

E^{f) ^ E^[Q) - x^f . (3) 

The two fit parameters used in this approach are the es- 
cape rate at zero force, ka — uj{0)e^^^''^'^/^'^'^ ^ and the 
dissociation length, Xh- Within Bell's model, it is impos- 
sible to uncouple the intrinsic time and energy scales of 
the system, because multiplication of w(0) by an arbi- 
trary constant and addition of k^T times the logarithm 
of that constant to £'b(0) leaves the force-free rate value, 
/cq, and hence the statistics of escape events, the same. 

One of the central problems of the present work is gain- 
ing more information about the system than Bell's model 
allows. Obviously, in order to make it possible, one needs 
to use a functional form of the escape rate k{f) involv- 
ing more fit parameters than in Bell's Ansatz. Indeed, 
it has been suggested [1, [3, H, [25, 26, 30.] that one can 
actually estimate the value of the force-free activation 
barrier i?b(0) if one makes a more realistic assumption 
about how the barrier decreases from this value with in- 
creasing force; the linear decrease ([3]) is replaced with a 
non-linear potential. Then, the force-dependent rate in- 
volves more than two fit parameters, the force-free barrier 
height £'b (0) being one of them. 

It is intuitively clear that the more rate parameters 
one has, the better one can characterize the system of 
interest. On the other hand, if the number of fit pa- 
rameters is too large, then not all of them may be de- 
termined sufficiently accurately from the experimental 
data. Therefore, a relevant question is: what is the max- 
imal number of rate parameters, which one can establish 
reliably? For a quite common and general parametriza- 
tion of the potential we show that the highest number 
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of model parameters one can determine from fitting the 
experimental data is three, while inclusion of additional 
parameters into the theory does not improve the quality 
of the fit. These three parameters are related to the force- 
free escape rate, the height of the activation barrier, and 
the dissociation length. Furthermore, we show that the 
resulting fit values of these parameters strongly depend 
on the assumption concerning the manner in which the 
activation barrier decreases with the applied force. In 
other words, one cannot determine them uniquely with- 
out having this information. 

II. THEORETICAL MODELING OF BOND 
RUPTURE 

For a quantitative analysis of dynamic force spec- 
troscopy experiments one usually assumes that the force 
f{t) increases linearly, 

fit) = KVt , (4) 

where v is the (constant) pulling velocity, and k, the rel- 
evant total elasticity of cantilever, linker molecules, re- 
ceptor and ligand, cf. Fig. 1. Using the fact that the 
extension s of the complex at time < is s = we see 
that ^ is equivalent to the assumption of a linear force- 
extension characteristics. 

Combining Eqs. and we arrive at the central 
experimental quantity, namely the probability density of 
escape events at the force value / for a given velocity v 
and a given set of parameters fj. characterizing the escape 
rate k{f): 

MfM^-'-^^'-Mle-^^fy-, (5) 
aj Kv 

g{I):^ f dfkin ■ (6) 

For later convenience, the rupture force distribution ^ 
is written as a conditional probability, conditioned on the 
values of the model parameters n and the pulling velocity 
V. While, in practice, the latter is directly accessible from 
the measurement, the model parameters jj, have to be 
inferred from the distribution of rupture forces. 

It is possible to evaluate the integral ^ numerically 
for any functional form of the rate k{f ). However, the 
exponentially increasing character of this function, see 
Eq. ^ , allows one to derive a very efficient analytical ap- 
proximation for the integral. Namely, in view of Eq. ([1]), 
the main contribution to the integral g{f) comes from 
the /'-region just below /. This allows one to expand 
lnfc(/') in the vicinity of the upper limit of integration / 
to the first order [T7| : 

fc(/')-fc(/)e^^(^)(^'-^) , (7) 
A„(/):=(-ir-^^^. (8) 



With this approximation, the integral in ([6]) is given by 
ff(/)-ffi(/) = |^(l-e-^^(^)^) . (9) 

This simple approximate formula may be sufficient for 
most practical purposes. However, if the deviations of 
lnfc(/') from linearity are important, one can use the 
second-order approximation for this function near /' = /: 

fc(/') ^ fc(/) e^i(/)(/'-/)-A2(/)(/'-/)V2 ^ (10) 
allowing one to evaluate the rate integral in ^ as 

9{f) « 9.{f) = (11) 

where erf (a;) := 2tt^^^'^ dy e^^ is the error function. 
Additional numerical analysis has shown that the inac- 
curacy of this expression is smaller than 1% for all rea- 
sonable functional choices fc(/) which we have checked. 

III. PARAMETER ESTIMATION: PROPERTIES 
OF THE MAXIMUM LIKELIHOOD ESTIMATOR 

Let us assume that a specific model and thus the par- 
ticular form of the probability density pi in ([5]) can be 
considered as given. Then the remaining task is to esti- 
mate the model parameters /x from a given set of A'^ rup- 
ture forces f = {fi}iLi and pulhng velocities v = {vi}^^. 
There exist different "recipes" for doing this, called es- 
timators. Each of them can be formally represented by 
some function /i(f), indicating that the estimate, being a 
function of the random variables f , is a random variable 
itself. Now, the main question is: what is the optimal 
estimate of the model parameters that can be extracted 
from the given set of N rupture forces? Stated differ- 
ently: which recipe yields estimates /i(f ) of the model 
parameters which are on average over many data sets f 
closest to the "true " model parameters! 

In this section we discuss some properties of estimators 
and show that under realistic experimental conditions, 
given in single-molecule pulling experiments, no estima- 
tor /i(f ) outperforms the maximum likelihood estimator. 
The reader who is not interested in the mathematical 
details may skip the subsequent subsections A-C and im- 
mediately proceed to subsection D, where we summarize 
the main steps which are necessary for a practical appli- 
cation. 



A. Maximum Likelihood Estimator 

Our starting point is the probability p to observe a 
given set of A^ rupture forces f = measured at 
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pulling velocities v = {vi}fLi. Since the fi are statisti- 
cally independent, this probability reads: 

N 

pii\(i,^)=l[pi{f,\fM,v,) . (12) 

i=l 

The proceeding consists in simply maximizing (jl2p with 
respect to /i [23, HO] ; usually, this has to be done numer- 
ically. For any given f and v the corresponding set of 
parameters fi* = fi*(i,v) is called the maximum likeli- 
hood estimate. 

Intuitively, the properties of this estimator are most 
easily understood within the framework of Bayesian infer- 
ence [31I, [S^l ■ The quantity in is called likelihood and 
plays a central role in the Bayesian approach. Extend- 
ing the notion of "probability" in the sense of "degree of 
belief" to the model parameters /x, the joint probability 
p(f , fi, v) can be written in terms of conditional proba- 
bilities p{...\...) either in the form p(/i.|f, v)p(f, v) or in 
the form p(f l/x, v)p(/x, v), yielding Bayes' theorem: 

p(/x|f,v) =p(f|/x,v)p(/x,v) b(f,v)]-i . (13) 

The left hand side represents the "likeliness" of fi, given 
the data f , v, and hence is clearly of central interest for 
our purposes. Considering also the right hand side as a 
function of /x, it is equal to the likelihood from (fT^ times 
the so called prior probability p(fi,-v), encapsulating all 
our knowledge about /x before the measurement, times 
a /x-independent factor. Thus, determining p(/x|f, v) by 
means of the right hand side of (fT5|) provides the central 
"recipe of learning" within the Bayesian approach (sil . 

Regarding actual practical application of Bayesian in- 
ference, the determination of the prior probability is the 
most problematic point. Different recipes for selecting 
an appropriate prior exist. Common choices are distri- 
butions which are uniform in the parameters or the log- 
arithms of the parameters. Rigorous justifications are in 
general not possible and one is left with postulating some 
heuristic ad hoc Ansatz. 

However, dynamic force spectroscopy usually provides 
large data sets, i.e. large TV. Then the likelihood (fT^ 
develops a narrow peak in the region of its maximum fi* 
(see next section) and the prior p{fj,,'v) in (jl3p . though 
usually unknown in detail, can be considered as approx- 
imately constant, i.e. p{fi\i,v) c>c p(f |/x, v). Given f and 
V, the likelihood p2)) thus quantifies the "likeliness" that 
the "true" model parameters are /x. 

The upshot of the above intuitive considerations is that 
maximizing (jl2p with respect to /x should yield the best 
possible guess for the unknown true parameters. Further- 
more, the statistical uncertainties of this estimate should 
be somehow related to the width of the likelihood. In the 
following subsection, we leave this intuitive level and turn 
to a more rigorous discussion of the asymptotic proper- 
ties of the maximum likelihood estimator. 



B. Asymptotic Properties 

Let us assume that the rupture forces fi have been 
sampled according to the "true" distribution pi (/^ | /Xq, v^) 
with unknown, "true" model parameters /Xg. For a given 
set of rupture forces f and pulling velocities v the max- 
imum likelihood estimate can then be determined as de- 
scribed above. Upon repeating the entire set of iV pulling 
experiments with the same set of pulling velocities v, a 
different set of rupture data f will be sampled, yielding 
a different maximum likelihood estimate /x*. While the 
probability distribution of f is given by (fT^ with fi — (Iq, 
what can we say about the distribution of the maximum 
likelihood estimates /x*? 

To answer this question we first exploit the fact that 
in typical single-molecule pulling experiments for each 
pulling velocity several hundred rupture forces are mea- 
sured. The resulting set of rupture force data f is thus 
quite large and it is convenient to rewrite the likelihood 
(HH) as 

p(f|/x,v) = exp{-iVsw(f,/x,v)} (14) 

N 

sjv(f,M,v) := -iV-i^lnpi(/,|/x,w,) . (15) 

i=l 

Furthermore, we assume that rupture forces have been 
measured at Z different pulling velocities vp, (3=1, Z , 
and that the relative frequency with which the different 
pulling velocities v are sampled, converges towards a well 
defined limit p{v) = J2p=i Pf^H'^ — Vfs) for N ^ 00. Then 
it follows from the law of large numbers [s^ that 

SAr(f,/x,v) s(/x) := -(lnpi(/|/x,'y))i (16) 

for N —^ 00, where (• • •)! indicates an average over / and 
V with weight pi(/|/Xq, v) p{v). Hence, sat is an intensive, 
entropy-like quantity. Observing that s(/x) — s(/Xq) is a 
relative entropy of the form (ln(pi(/|/XQ, u)/pi(/|/x,, ti))i, 
and using the fact that pi(/|/x, v) is normalized with re- 
spect to / for every choice of the parameters /x, we obtain 

s(/x) ~ s(/Xo) = (17) 

Jdvp{v) J dfpiif\fj,,v)[RlnR-R+l] , 

with R := p(/|/Xq, u)/p(/|/x, u). Finally, using the in- 
equality 

0<J dx \nx = R\nR- R+ 1 Vi? , (18) 

we see that s(/x) — s(/Xq) > and that s(/x) — s(/Xq) = 
if the expression in the square brackets on the right hand 
side of (|17p vanishes for all /. Thus, s(/x) has a unique 
absolute minimum at /x = /Xq [33|. Since sn converges 
for large N toward s according to (|16[) , also the minimum 
/X* of the former converges to the minimum /Xq of the 
latter, i.e. the maximum likelihood estimate is a so-called 
consistent estimate 1341. 
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For large, but finite N values, /i will be close to /i.*. 
Consequently, we can expand SAr(f,/i.,v) up to second 
order about its minimum at fi* and neglecting terms of 
order 0{1/\/N), the Hessian matrix of s n {{ , fj,* , -v) can 
be replaced by the Hessian H — H{^q) of s(/^o)' which 
is generically positive definite, i.e. 

sw(f,M* + A,v) = SAr(f,/x*,v) + Ati?A/2 . (19) 

For large N this is a very good approximation for all /x- 
values, and p(f|/i.,v) approaches a very sharply peaked 
Gaussian about /^*, 

p(f|/x,v) cxexp{-A^(/x-/x*)ti7(/x-/x*)/2} . (20) 

We can now determine the first moments of the dis- 
tribution of the maximum likelihood estimate fj.* (upon 
many repetitions of the same experiment). Differentiat- 
ing and choosing A ^ fiQ — fi* , results in 

/x*-/Xo = -i7-iasjv(f,/Xo,v)/a/x . (21) 

Averaging over f yields zero in the right hand side, as 
can be inferred from p^ . and the fact that fi^ is 
the minimum of s. Hence, 

(m*> = Mo , (22) 

where (• • •) indicates an average over f with weight 
p(f I/Xq, v) for a given set of pulling velocities v. Equation 
(1^^ thus shows that the maximum likelihood estimate is 
"unbiased" for large N. 

With (PT|) the determination of the second moments is 
straightforward. Using 

^ A,^(f,/,^,v) A,^(f,/,„,v)^ ^ ^^^^ ' (23) 

gives the covariance matrix for the maximum likelihood 
estimate: 

([/x*-/x„][/x*-/Xo]t> = (7ViJ)-i . (24) 

Observing that {N H)~^ is also the covariance matrix of 
the distribution from ^U\j we arrive at our 

First main conclusion: For any given, sufRciently large 
data set f, the expected deviation of the concomitant 
maximum likelihood estimate fi* from the "true" param- 
eters fiQ immediately follows from the "peak-width" of 
likelihood p^ . considered as a function of fi. 

Similarly, using the central limit theorem, one can 
show (see Appendix . that fj,* is Gaussian distributed, 
yielding with (pO|) our 

Second main conclusion: Apart from the peak position 
and a normalization factor, the likelihood for one 
given data set f looks practically the same as the dis- 
tribution of the maximum likelihood estimates fi* from 
many repetitions of the N pulling experiments. 



C. Cramer-Rao bound 

It should be noted that, in order to derive the above 
two main conclusions, we did not make any use of the 
Bayesian formalism P^ at all. The latter only served to 
acquire an intuitive idea about the meaning of the likeli- 
hood p^ . At this intuitive level, we have seen that the 
left hand side of p^ is very well approximated by the 
sharply peaked Gaussian in and hence it is reason- 
able to expect that its maximum fi* should be the best 
possible guess for the unknown true parameters /Xg that 
possibly can be inferred from a given set of data f . A 
more rigorous line of reasoning starts with an arbitrary 
"recipe" /i(f ) of estimating the true parameters /Xq from 
a given data set f. The only assumption is that this 
recipe is unbiased, i.e. upon repeating the same exper- 
iment many times, the "true" parameters are recovered 
on average, (/j(f)) = P-n- By generalizing the well-known 
Cramer- Rao inequality [33|, which in turn is basically 
a descendant of the Cauchy-Schwarz inequality, one can 
show [sS] for any such "recipe" /i(f) that 

([/i-Mo][/i-/Xo]^) -(A^i/)-^ >0 , (25) 

i.e. the matrix in the left hand side is non-negative defi- 
nite. Comparison with yields our 

Third main conclusion: There is no unbiased estima- 
tor fi of the true parameters /Xg which on the average 
outperforms the maximum likelihood estimate fi*. 

The remaining possibility that a biased estimator may 
be even better is rather subtle to treat rigorously, but 
intuitively this seems quite unlikely. Furthermore, in the 
above conclusion we exploited the relation ([M]) which is 
strictly correct only for asymptotically large N. Finally, 
also the criterion of minimizing the left hand side in (j25p 
itself is in principle debatable, but hardly in practice. 



D. Parameter Inference: Main steps for the 
practical application 

We now briefly summarize the main steps of the max- 
imum likelihood method for evaluating single-molecule 
pulling experiments. The first step consists in specifying 
the dependence of the rupture force probability pi(/|/x, v) 
on the model parameters /x and the pulling velocity v 
within any given theoretical description. Then for the set 
of rupture forces f — {fi}fLi f^nd corresponding pulling 
velocities v = {vi}fLi the logarithm of the likelihood pi)) 

N 

-Ars^(f,At,v) = ^lnpi(/,|/x,t-0 (26) 

1=1 

is maximized with respect to the model parameters fi. 
Usually this step has do be accomplished numerically. 
The position of the maximum defines the most probable 
parameters fi* which are on average closer to the true 
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model parameters than any other estimate. The statisti- 
cal uncertainties of the parameters can then be estimated 



(27) 



where Hn denotes the Hessian matrix of sat evaluated at 
the most probable parameters /x* and /Xq the true model 
parameters. Moreover, the distribution of the maximum 
likelihood estimate is Gaussian with mean /Xq. 



IV. APPLICATION TO SINGLE-MOLECULE 
FORCE SPECTROSCOPY: BELL'S MODEL 

Combining the result ([5]) with approximations ([7]), ^ 
[or (fTO|) . (|lip ]. one can apply the maximum likelihood 
approach from the previous section to deduce the rate 
parameters for any exponentially increasing escape rate 
(fT|). In what follows, however, we will focus on two effi- 
cient rate approximations [see Eqs. ((28|) and ([35| below], 
which allow one to evaluate the integral from Eq. ([6|) 
analytically. 

The first common approximation, originally due to Bell 
[Tl| . consist in the linearization of the force dependent 
potential barrier according to ^ and in neglecting the 
force dependence of the pre-exponential factor uj{f) in 
the Kramers rate ([Ij, resulting in [cf. ([3]) and Fig. 2] 



k{f) = fcoexp 



fcRT 



exp(A + af) 



(28) 



where ko := (jj(0) exp(— i?b(0)/fcBT) is the force-free dis- 
sociation rate, Xb the dissociation length (distance be- 
tween potential well and barrier), and A :~ Infco: Q! := 
Xh/k^T are convenient abbreviations. 

Substituting Eq. ((28|) into (O, ^ and going over to 
/ as independent variable, a straightforward calculation 
yields the probability density of rupture events for Bell's 
model: 



■ exp 



(29) 



The rupture force density (j29p is conditioned on fi — 
(A, a), and v. As usual, we assume that the pulling ve- 
locity V is known exactly for each measurement, and sim- 
ilarly for the elasticity k appearing on the right hand side 
of (P^ . The remaining model parameters to be estimated 
from a given set of rupture forces f = {fi}fLi measured 
at pulling velocities {vp}^^^ with relative frequencies pp 
are therefore fj. = (A, a). 



A. Statistical uncertainties of Bell's model 

For the above specified model we can calculate s(/x) = 
— (lnpi(/|/x, i;))i as defined in (|16p analytically, if the 
dimensionless quantity Tp := e'^° / {nvpao) is small for 



all pulling velocities. In fact, < 1 is equivalent to 
the assumption that the distribution of rupture events 
Pi(/|/ig, i;/?) has a maximum at some force /| > 0. De- 
tails of the calculations are given in Appendix [Bl result- 



ing m: 



s{^l) 



z 

-\ + Y,Ppsp{ti) + 0{Tp) , (30) 

(3=1 



S/3(m) := ln(KW/5) -I- ri{C + Inr^j) + 



F (r; -|- 1) 
KVpa [rpY' 



with ?] := a/ao, C w 0.577 Euler's constant, and F(-) the 
Gamma function. Differentiating ([50]) twice with respect 
to the model parameters, a straightforward calculation 
yields the Hessian H — H{fi(f) of s(/1q). Finally, invert- 
ing this (2 X 2) matrix, we obtain the variance of the 
maximum likelihood estimate of the two parameters: 



{ia*-aor) 



((A*-Ao)2) 



an 



1 13=1 ^ 

N f +a2(ln(r)) 



(31) 



(32) 



where 



z 

a2(ln(r)) = ^p^W(r^) 

f3=l 



Z 

(3=1 



(33) 



is the variance of the logarithm of the loading rate r := kv 
and 



ifp) = -l/ao{C + \n{e^° /{rpao))) 



(34) 



the expected rupture force at loading rate rp (again ne- 
glecting terms of order 0{Tp)). For practical application 
of Eqs. ([?T|) and ([5^ the true model parameters /Xq in the 
right-hand sides of the equations have to be replaced by 
the concomitant maximum likelihood estimate fi* . For 
large N this is a very good approximation. 

In single-molecule pulling experiments only a limited 
range of pulling velocities v is accessible, i.e. vp € 
[i^miiD 'I'max]- Now the question arises: for which distri- 
bution of pulling velocities p{v) do the statistical uncer- 
tainties (|3ip . ((32)) of the estimated parameters become 
minimal? Recognizing that the variance of the maxi- 
mum likelihood estimate of a depends on the distribu- 
tion of the pulling velocities solely via the term (T^(ln(r)) 
in the denominator of (j31[) . while this distribution en- 
ters into the expression for the statistical uncertainties 
of A* also via the terms {fp) in the numerator of ([5^ . we 
see, that it is, in general, not possible to simultaneously 
minimize the two uncertainties. Given a fixed number 

of pulling experiments, the minimization of the vari- 
ance {{a* — ofg)^) of a* is equivalent to a maximization of 
(T"(ln(r)). Under the constraint vp € [v 



mull ^maxj 



for all 
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pulling velocities V[}, this maximum is obviously reached 
if half of the rupture forces have been sampled at a pulling 
velocity as large as possible, i.e. Wmax, and the other half 
at a pulling velocity as small as possible, i.e. Wmin- This 
result is independent of the values of the true model pa- 
rameters fiQ. Regarding the distribution of A*, the situa- 
tion is more complicated. Given the analytic expression 
(15^ for the error, the calculation is straightforward, but 
not shown here for the following two reasons. The "best" 
choice consists again in sampling just at the two extreme 
pulling velocities v^in and Wmax- The relative number of 
pulling experiments for each of the two pulling velocities 
is, however, non-trivial and depends on the values of the 
true model parameters /Xq. It is, therefore (in contrast 
to the result for a*), only of limited use for a real exper- 
iment. 



B. Illustration for computer generated data 

We now illustrate the findings from sections IIII Bl and 
IIII CI To avoid uncontrollable experimental inaccuracies 
and uncertainties regarding the "true" model and the 
"true" model parameters /Xq, we numerically generated 
synthetic rupture data f by "simulating" an idealized ex- 
periment on the computer according to the probabilistic 
"laws" (|28p . with given parameters fj, — fj,Q. Hence 
all remaining uncertainties are statistical finite N effects. 
Fig. [3] shows the results for 10000 repetitions of a com- 
puter experiment, each sampling N = 400 rupture forces 
f according to (j29p with experimentally realistic param- 
eter values ao = 0.1 pN~^ and Aq — —5. Since two- 
dimensional distributions are difficult to compare graph- 
ically, we focus on the marginal distributions. For each 
of the 10000 experiments, the maximum /x* = (A*, a*) 
of the likelihood p2|) . considered as a function of fi, was 
determined numerically. The distribution of the result- 
ing A*- and a*-values are depicted as histograms in Fig. 
[31 The standard deviations of the maximum likelihood 
estimate for A and a, determined from the 10000 experi- 
ments are s(A*) = 0.20 and s{a*) = 0.0023 pN"! respec- 
tively. These values coincide with those obtained from 
the analytical approximations pT|) - (|34p up to the third 
non- vanishing digit. Replacing the parameters /Xq on the 
right-hand side of (pij) . (|32p by the maximum likelihood 
estimate fi* for one given data set thus provides reliable 
estimates for the statistical uncertainties. 

Furthermore, for the first 15 of the 10000 experiments, 
after integrating over the other parameter in the likeli- 
hood (fT^ , shifting the peak position from /x* to iJi^o , and 
normalizing, the resulting marginal distributions were 
plotted in Fig. [31 They closely agree with the histograms. 
These observations illustrate very convincingly our two 
conclusions from section llll Bl above. In particular, finite- 
N corrections are apparently very small for the typical 
parameter values used in this example. Also the prac- 
tically perfect Gaussianity of the distributions is as ex- 
pected, cf. (HOI)- 



(a) 



(b) 




0.09 0.1 0.11 

a' [pN"'] 

FIG. 3: Solid histogram: Distributions of the first and second 
components of the maxima /x* — (A*, a*) of the likelihood 
(|12p for 10000 "computer experiments". For each of them, 
= 400 rupture forces / were sampled according to ((29}, 
100 for each of the 4 loading rates kv = 50, 200, 1000, 5000 
pN/s and with "true" parameters Ao = —5 and ao ~ 0.1 
pN"\ These are typical numbers in "real experiments" [J]. 
For sake of better visibility the bin-width of the histograms 
is much larger than the optimal bin width for a Gaussian dis- 
tribution (Appendix [Cj . Thin lines: Likelihood p2|l for the 
first 15 of the 10000 experiments after integrating over the 
other component, shifting the maximum to f^i,o, and normal- 
izing (some are almost indistinguishable). Dotted histogram: 
Distribution of the estimates for A and a according to the 
"standard method", as described in the main text. 



Let us finally compare the performance of the maxi- 
mum likelihood estimate with that of the most widely 
used "recipe" of parameter estimation in the field of 
single-molecule pulling experiments. This consist of the 
following steps: (i) Fit a Gaussian to the observed rup- 
ture force distribution for a fixed pulling velocity v and 
approximate the most probable rupture force /* by the 
maximum of that Gaussian, (ii) Plot /* for different 
V versus ln(i;) and fit the resulting points by a straight 
line, (iii) Assume that the model ([^S)) . (|^^ is applica- 
ble and deduce its model parameters fi — (A, a) from 
the slope and the axis intercept of the straight line as 
detailed e.g. in P, [H, [H, [H, [lH . We have applied this 
procedure to each of the 10000 experiments in Fig. [31 and 
plotted the distribution of the resulting estimates for A 
and a in Fig. [3l The systematic bias of the estimate 
for A can be traced back to fitting a Gaussian, which is 
symmetric about its maximum, to an asymmetric "true" 
distribution (|29p [ssj . while the suboptimal variance of 
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the estimate for both A and a signals that quite some in- 
formation is fost by only considering the most probable 
rupture forces /*. Hence, the maximum likelihood esti- 
mate represents a substantial improvement compared to 
the so far "standard method" of data evaluation in this 
field. This is in agreement with our conclusion from sec- 
tion llll Cl We have also directly compared the maximum 
likelihood estimate with other known "recipes" of evalu- 
ating single- molecule rupture data, e.g. [S^. In all cases 
we found that the maximum likelihood was superior. 

V. EXTENSION OF BELL'S MODEL 
A. Rate Ansatz 

As can be seen from Fig. [21 a linearization (128]) of the 
force dependent potential barrier Eh{f) is relatively good 
for small forces. For a larger forces, the distance between 
potential extrema decreases, leading to a weaker sensitiv- 
ity of the barrier height to force variations upon further 
puUing than in Eq. (^5]) . 

Models including this effect, in general, rely on some 
assumptions concerning the shape of the energy land- 
scape. Typical choices are Morse potentials, harmonic 
potentials with a cusp barrier, two parabolas at the po- 
tential extrerna joined at a rnidpoint, and linear-cubic 
potentials [1, M [11, [13, S H S lip ■ It has been 
suggested in j25| that for sufficiently high barriers, i.e. 
forces substantially smaller than the critical force, the 
dissociation rate can be written in a unique form 

k{f) = (1 - 7a//e)i/''-i e^+'^[i-(i-'^"^/<^)'^^l (35) 

with three model parameters fi = (A, a, e) and fixed ex- 
ponent 1/7. Here, A and a have the same physical mean- 
ing as in Eq. (^5)) . and e := Eb{0)/kBT stands for the 
force- free activation energy barrier in units of the thermal 
energy fceT. 

The extra parameter 7 controls the manner in which 
the barrier height decreases with the applied force. We 
note that, physically, this parameter should be in the 
range 76 (0,1] since 7 < would imply a positive barrier 
for all / > 0. On the other hand, the first derivative of 
the barrier height equals minus the distance between the 
potential extrema corresponding to a given force value. 
Since we expect this distance to decrease with the force, 
we conclude that the second force derivative of the barrier 
height must be positive, excluding 7-values greater than 
1. Specifically, for 7 = 1 the parameter e drops out and 
one recovers Bell's model (|28p . 7 = 2/3 reproduces the 
Kramers rate for a cubic reaction potential, and 7 = 
1/2 corresponds to a parabolic potential well with a cusp 
barrier. 

Substituting Eqs. (g]), ([23) into (O, one derives the 
survival probability of the bond up to force /: 

n(/)=exp - — ^ , (36) 

\ Kv a I 



and the probability density of rupture events follows from 
Eq. ©. 

As already mentioned, the application of Kramers re- 
action rate theory requires that the potential barrier 
£'b(/) be sufficiently high (compared to the thermal en- 
ergy k-QT). Thus, all approximations are only valid 
for forces substantially smaller than the critical force at 
which the barrier vanishes, /c = 6/(70;). 

The above discussion suggests that dynamic force spec- 
troscopy should, in principle, provide the possibility, 
not only to determine the force-free dissociation rate 
/co = exp(A) and the dissociation length a;b, but also 
the force-free activation energy barrier i?b(0) = tk^T. 
Naturally, the question arises how accurate these esti- 
mates will be and whether the inferred values critically 
depend on the chosen theoretical model, in particular on 
the value of the parameter 7. These questions will be 
addressed next. 



B. Numerical experiment 

There is an ongoing debate in literature about which 
of the three exponents, i.e. 7 = 1, 7 = 2/3, or 7 = 1/2, 
is most appropriate to use when evaluating experimental 
rupture data SI, [H [H, HI] . Taking for granted that 
one of the three models approximates the "truth" satis- 
factorily, choosing fi ^ fi* is ~ according to our above 
conclusions - the closest one can get to the "full truth" 
on the basis of one given data set f . In case of disagree- 
ment about the "true" 7-value, a fully objective selection 
criterion seems unavailable in principle. In practice, the 
usual criterion is the comparison with the basic "true" 
quantity observed experimentally, namely the distribu- 
tion of rupture forces. 

For the example shown in figure U] we sampled N ~ 
2000 rupture forces f according to the distribution ([5]) 
with rate psp . In the simulations, we have taken 7 = 2/3 
and realistic parameter values Aq = — 5, ao = 0.1 pN"-'^, 
eo = 15. 

In order to test the maximum likelihood method, for 
the resulting data set we determined the estimate fi* for 
several possible 7-values. Since in a real experiment, the 
value of the exponent 7 is not known a priori, also during 
the fitting, the 7-values used were not necessarily coin- 
cident with the "true" value used for data generation. 
The bin- width of the histograms in Fig. |4] was chosen as 

— 1/3 

hf3 = 2.83sNpN„ , where sat^ is the standard deviation 
of the rupture forces measured at pulling velocity vp . As 
discussed in Appendix[Cl this is the optimal choice of the 
bin-width for the Bell model; although this model is in 
fact not the true one in our numerical experiment, this 
choice of the bin-width remains suitable, as every distri- 
bution separately can be very well fitted with Eq. (|29)) . 

Comparing in Fig. |4] the resulting distributions 
Pi{f\fj,*,v) for three different 7 values, we observe the 
following. Already the two-parametric Bell's Ansatz 
(7 = 1) reproduces the experimental distribution of rup- 
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FIG. 4: (Color online) Rupture force distribution for different 
loading rates kv. Histograms: numerically generated rupture 
forces according to ((5]), ((35|) with 7 = 2/3, Ao = —5, ao = 0.1 
pN~^, eo = 15. For each kv, we sampled 500 forces, i.e. = 
2000. The bin- width is chosen according to eq. (|C6[l . The 
maximum likelihood fits according to (|35p for 

7 = 1/2 and 7 = 2/3 (red solid) are not distinguishable within 
the line width in this plot. For 7 = 1/2, the fit parameters 
have the following values: A = —5.37, a = 0.110, and e = 
17.2. For 7 = 2/3, the fit results are slightly closer to the 
true parameter values, namely, A = —5.20, a — 0.104, and e = 
14.64. Blue dashed distribution: same for Bell's Ansatz 7 = 1, 
with fit results A = —3.81 and a = 0.072. Upon repeating the 
entire "numerical experiment" , the resulting plots always look 
practically the same. 



ture forces with reasonable accuracy, see the dotted line 
in Fig. m However, if one increases the number of fit 
parameters to three, one obtains the fit distributions no- 
tably different from the Bell's curves, see the solid lines 
corresponding to 7 = 1/2 and 7 — 2/3. We also note that 
the difference between the curves corresponding to these 
two values of 7 is smaller than the line thickness. This 
means that if one treats the exponent 7 as a fourth fit pa- 
rameter, then its precise value cannot be determined by 
fitting the experimental rupture force distribution. With 
respect to the three remaining fit parameters, their val- 
ues are rather close to each other for the fits with 7=1/2 
and 7 — 2/3, with the largest discrepancy between the 
fit values of the force-free barrier height (see the caption 
in Fig.HD. 

When fitting the real experimental data, one does not 
know a priori the true value of the exponent 7. It is 
therefore of interest to find out, how the remaining fit 
parameter values depend on the assumption with respect 
to this quantity. Fig. [5] shows the fit results obtained 
for different assumed values of 7 within the physically 



meaningful range. All of the fitting curves obtained for 
different 7-values from Fig. [5] coincided within the line 
thickness. We observe that the resulting fit values of 
the force-free escape rate and the dissociation length 
k^Ta are not very sensitive to the choice of the exponent 
7. At the same time, the force-free barrier height value 
inferred from the fit depends approximately linearly on 
the choice of this parameter, and can assume values dif- 
fering by as much as a factor of 2 at extreme 7-values. 
This means that when fitting the experimental data, the 
value of the force-free barrier height will be determined 
with the least accuracy. 

In view of these observations, an interesting question 
arises: are the approximations (|24p and pT|) - (|34p for 
the statistical uncertainties of the model parameters still 
valid for those models, about which we (in our case) 
know that they are not true? To study this point, we 
have repeated the above described procedure for 10000 
data sets, each generated in the same way and with same 
"true" parameters as the data set shown in Fig. [4l The 
distributions of the inferred parameters fi* for the three 
different 7-values are depicted in Fig. [5] as histograms. 
Following section [Til Bl they should be bell-shaped with 
variance given by Eq. ([24]). We have evaluated this ex- 
pression at the mean values of the inferred parameters fi* 
for each of the three 7-values respectively. For 7 = 2/3 
and 7 = 1/2 this had to be done numerically, whereas for 
7 = 1, Eqs. dnil-llMl) could be employed. The result- 
ing distributions, centered about the mean value of the 
inferred parameters, are shown in Fig. [S] as solid lines. 
For A and a they closely agree to the histograms. With 
respect to the parameter e, the approximated variance 
agrees very well with the empirically determined, but 
finite-iV corrections to the full distributions which are 
not symmetric about their center are apparent. 

In conclusion, although the models for 7 = 1 and 
7 = 1/2 are not the true models, equation (|24p still 
yields very good approximations for the statistical un- 
certainties of the parameters, i.e. for their distribution 
upon repeating the same experiment many times. Never- 
theless, comparing these distributions to the distribution 
of the parameters for the "true" gamma- value 7 = 2/3, 
we see, that by choosing the wrong model, the system- 
atic deviations to the parameters are much larger than 
the statistical uncertainties. As one will, in general, be 
unsure about the true underlying energy landscape, and 
thus about the "true" 7-value, this point is essential if 
one wishes to use the inferred parameter values in an- 
other context than the interpretation of single-molecule 
pulling experiments. 

We would like to mention that the deviations between 
the distributions resulting from the different 7 values in- 
crease with the pulling velocity [2^. Hence, by increas- 
ing the range of accessible pulling velocities [I'min, I'max] a 
clearer distinction between the models is possible. How- 
ever, for precise measurements with the AFM the loading 
rate kv is limited to a few orders of magnitude, compa- 
rable with our values. 
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FIG. 5: Maximum likelihood fit values for the data set from Fig. |4] fitted with rupture force distribution ([Sjl supplemented 
with the escape rate (|36ll . Each data point was obtained by fitting the same data set, but assuming a different value of the 
parameter 7. 
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FIG. 6: (Color online) Histograms; Numerically determined distribution of the single components of the maxima fi* = 
(A*,a*,e*) of the likehhood (HIJ, ^ with 7 = 1/2 (green), 7 = 2/3 (red) and 7 = 1 (blue) for 10000 "computer experiments". 
For all data sets the rupture forces were generated numerically according to ([5]), (|35p with 7 — 2/3, Ao = —5, ao = 0.1 pN~^, 
£() = 15. For each kv, we sampled 500 forces, i.e. A*' = 2000. The bin-width of the histograms is see 
Appendix[Cj. Solid lines; Gaussian approximations to the distributions with covariance matrices ((24}. They have been shifted 
so that their maximum coincides with the mean of the empirically determined distributions. For sake of better visibility the 
distributions for 7 = 1 have been rescaled by an appropriate value. 



VI. CONCLUSIONS 

In this work we have shown that the maximum likeli- 
hood approach is an extremely simple, general, and pow- 
erful method for parameter estimation in the contect of 
single-molecule force spectroscopy. For large data sets it 
outperforms all other estimates. Furthermore, approxi- 
mations to the statistical uncertainties of the parameters 
are available once the parameters are estimated. In the 
case of the standard Bell model we were able to derive an 
analytical expression for these uncertainties in terms of 
the model parameters and the distribution of the applied 
loading rates. For more general models, the uncertainties 
can be determined numerically. 

When fitting the experimental data, one usually adopts 
some functional form of the force-dependent escape rate 
involving several fit parameters. By means of a numerical 
example, we have demonstrated that the largest number 
of such parameters that can be determined from the ex- 
periment is three. These parameters are related to the 
force-free value of the rate, the dissociation length, and 
the barrier height in the absence of the force. Further- 
more, when fitting the experimental rupture force dis- 



tributions, one needs to make an additional assumption 
about the manner in which the escape rate decreases with 
the applied force. While the fit values of the force-free es- 
cape rate and the dissociation length depend only weakly 
on this assumption, the value of the force-free barrier 
height can be determined much less reliably. We have 
shown that even if the model adopted for the descrip- 
tion of the experiment is not the true one but predicts 
distributions of rupture forces similar to the measured 
distribution, the statistical uncertainties found from the 
maximum likelihood method very well approximate the 
dispersion of the estimated parameters upon repeating 
the same experiment many times. Often these uncertain- 
ties are much smaller than the systematic error resulting 
from choosing the "wrong" model. 
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APPENDIX A: ASYMPTOTIC DISTRIBUTION 
OF THE MAXIMUM LIKELIHOOD ESTIMATE 

In this section we proof the second main conclusion 
from section IIII B) namely that the distribution of the 
maximum likelihood estimate n* is Gaussian and looks, 
apart from the peak position, the same as the likelihood 
(pO)l for one given data set f . 

To keep things as simple as possible, we assume that 
there exist (small) integers for each of the Z pulling 
velocities vp so that np/rij = pp/ p^- Denoting by n the 
sum of all n/3, the total number of rupture forces can be 
written as TV = N'n and the set of N rupture forces f can 
be divided into N' subsets f*^ = {/f^} where for fixed fc, 
ffjj is one out of the np rupture forces sampled at pulling 
velocity vp. 

Having introduced this notation, we define N' new ran- 
dom variables 

X''^ -X'=(f'=) = -^^i/-l9/a/xlnpl(/5l/Xo,^;^) . 

/3=1 j=l 

(Al) 

We know already from our discussion in section lTlI BI that 
these random variables have an expectation value zero 
and a covariance matrix 



n 



(A2) 



Then, from the central limit theorem |34l it follows that 



1 ^' 



(A3) 



fc=i 



is Gaussian distributed with mean zero and covariance 
matrix 

([/.* - Mo] [m* - /^o]^) = ^H-' = ^H-' ■ (A4) 



TV 



APPENDIX B: STATISTICAL UNCERTAINTIES 
FOR BELL'S MODEL 

In order to determine the covariance matrix of the 
maximum likelihood estimate for the basic model, in sec- 
tion |IVX] the quantity s{fi) as defined in (fT5| had to be 
calculated. Details of this calculation are given below. 

We first calculate 



E{a) 



dfe 



oo 

1 = , df e'^fpiifWv) (Bl) 





■ exp 



r ao 



for an arbitrary loading rate r = kv and a > — ao. It 
is convenient to rewrite equation (jBip using the dimen- 
sionless quantities t := e^°/(rQ!o) and = a/ao and to 



substitute t = t exp{aQf) yielding: 

/oo 
dt t"e 



(B2) 



Using that r ^ 1 in typical AFM pulling experiments, 
we derive at 



(B3) 



E{a) ^T-'^Tir, + l) + OiT), 



where r(-) denotes the Gamma function. Equation (|B3[) 
directly gives the expected rupture force: 

(/)i = 7^i?(«-0) = -— (C + In r) + 0(r) (B4) 



with C ~ 0.577 the Euler constant. Using equations 
IB4[) and the linearity of the expectation value we obtain 



-(lnpi(/|/x,r))i = 

-A + ln{Kv) + 7]{C + Inr) -f' 



1 

Kva t''i 



r(,? + i) (B5) 



Finally, equation (jBSp together with the definition of the 
quantity s(^t) yields the desired result ([50)) . 



APPENDIX C: OPTIMAL BIN- WIDTH FOR 
HISTOGRAMS 

Let pi be a probability density function with two con- 
tinuous and bounded derivatives. For a sample of size N 
the histogram estimate pi of pi is defined as 



Mf) 



Ajv(/) 



Nh 



N 



(CI) 



with AAr(/) the number of values falling into the bin of 
width Hn around /. Then one can show 36, 37] that for 
large sample sizes the integrated mean squared error 



IMSE 



df{{piif)-Piif)f)i 



(C2) 



considered as a function of the bin- width is minimized by 



h* 
ri]\ 



N 



Jdfip'AfW 



1/3 



^-1/3 



(C3) 



Following the same lines as in Appendix [B] we obtain for 
the Bell model: 



df{p[{f\^,v)r^la'(l + o(—)) . (C4) 
8 \ \Kva J J 



2e^ 
Kva ^ 

Inserting (jC4p into (jC3p yields the optimal bin-width 

1 



(C5) 



It should be noted that in the limit -^S— <C 1, the optimal 
bin-width depends solely on the sample size N and on 
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the parameter a which determines the width of the dis- 
tribution, but neither on the force-free dissociation rate 
/co = exp(A) nor on the pulhng velocity v. 

If one wishes to determine the optimal bin-width prior 
to parameter estimation, one may make use of the well 
known result ((/ — ~ tt'^ / {Qa^) and choose 

Hn = 2.83sAr7V^i/3 , (C6) 



where s n denotes the standard deviation of the measured 
rupture forces. This estimate is pretty close to the op- 
timal bin-width hj^.causs ~ ZAQsnN"^/^ of a Gaussian 
distribution [361] . 
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